library(sf)
library(dplyr)
library(ggplot2)
library(scales)
library(ggmap)
library(leaflet)
ak_regions <- read_sf("data/shapefiles/ak_regions_simp.shp")
st_crs(ak_regions)
## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
class(ak_regions)
## [1] "sf" "tbl_df" "tbl" "data.frame"
#plot(ak_regions)
Use different projection with Alaska Albers
ak_regions_3338 <- ak_regions %>%
st_transform(crs = 3338)
st_crs(ak_regions_3338)
## Coordinate Reference System:
## EPSG: 3338
## proj4string: "+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
plot(ak_regions_3338)
summary(ak_regions_3338)
## region_id region mgmt_area geometry
## Min. : 1 Length:13 Min. :1 MULTIPOLYGON :13
## 1st Qu.: 4 Class :character 1st Qu.:2 epsg:3338 : 0
## Median : 7 Mode :character Median :3 +proj=aea ...: 0
## Mean : 7 Mean :3
## 3rd Qu.:10 3rd Qu.:4
## Max. :13 Max. :4
pop <- read.csv("data/shapefiles/alaska_population.csv", stringsAsFactors = FALSE)
head(pop)
## year city lat lng population
## 1 2015 Adak 51.88000 -176.6581 122
## 2 2015 Akhiok 56.94556 -154.1703 84
## 3 2015 Akiachak 60.90944 -161.4314 562
## 4 2015 Akiak 60.91222 -161.2139 399
## 5 2015 Akutan 54.13556 -165.7731 899
## 6 2015 Alakanuk 62.68889 -164.6153 777
Coerce the sf object using st_as_sf
pop_4326 <- st_as_sf(pop, coords = c("lng", "lat"), crs = 4326, remove = F)
st_crs(pop_4326)
## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
head(pop_4326)
## Simple feature collection with 6 features and 5 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -176.6581 ymin: 51.88 xmax: -154.1703 ymax: 62.68889
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## year city lat lng population geometry
## 1 2015 Adak 51.88000 -176.6581 122 POINT (-176.6581 51.88)
## 2 2015 Akhiok 56.94556 -154.1703 84 POINT (-154.1703 56.94556)
## 3 2015 Akiachak 60.90944 -161.4314 562 POINT (-161.4314 60.90944)
## 4 2015 Akiak 60.91222 -161.2139 399 POINT (-161.2139 60.91222)
## 5 2015 Akutan 54.13556 -165.7731 899 POINT (-165.7731 54.13556)
## 6 2015 Alakanuk 62.68889 -164.6153 777 POINT (-164.6153 62.68889)
Change coord reference system
pop_3338 <- pop_4326 %>%
st_transform(crs = 3338)
Do the spatial join
pop_joined <- st_join(pop_3338, ak_regions_3338, join = st_within)
head(pop_joined)
## Simple feature collection with 6 features and 8 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -1537925 ymin: 472627.8 xmax: -10340.71 ymax: 1456223
## epsg (SRID): 3338
## proj4string: +proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## year city lat lng population region_id region
## 1 2015 Adak 51.88000 -176.6581 122 1 Aleutian Islands
## 2 2015 Akhiok 56.94556 -154.1703 84 6 Kodiak
## 3 2015 Akiachak 60.90944 -161.4314 562 8 Kuskokwim
## 4 2015 Akiak 60.91222 -161.2139 399 8 Kuskokwim
## 5 2015 Akutan 54.13556 -165.7731 899 1 Aleutian Islands
## 6 2015 Alakanuk 62.68889 -164.6153 777 13 Yukon
## mgmt_area geometry
## 1 3 POINT (-1537925 472627.8)
## 2 3 POINT (-10340.71 770998.4)
## 3 4 POINT (-400885.5 1236460)
## 4 4 POINT (-389165.7 1235475)
## 5 3 POINT (-766425.7 526057.8)
## 6 4 POINT (-539724.9 1456223)
pop_region <- pop_joined %>%
as.data.frame() %>%
group_by(region) %>%
summarise(total_pop = sum(population))
head(pop_region)
## # A tibble: 6 x 2
## region total_pop
## <chr> <int>
## 1 Aleutian Islands 8840
## 2 Arctic 8419
## 3 Bristol Bay 6947
## 4 Chignik 311
## 5 Cook Inlet 408254
## 6 Copper River 2294
pop_region_3338 <- left_join(ak_regions_3338, pop_region)
## Joining, by = "region"
head(pop_region_3338)
## Simple feature collection with 6 features and 4 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -2175328 ymin: 405653.9 xmax: 773820 ymax: 2383770
## epsg (SRID): 3338
## proj4string: +proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## # A tibble: 6 x 5
## region_id region mgmt_area geometry total_pop
## <int> <chr> <dbl> <MULTIPOLYGON [m]> <int>
## 1 1 Aleutia… 3 (((-1156666 420855.1, -1159837 41… 8840
## 2 2 Arctic 4 (((571289.9 2143072, 569941.5 214… 8419
## 3 3 Bristol… 3 (((-339688.6 973904.9, -339302 97… 6947
## 4 4 Chignik 3 (((-114381.9 649966.8, -112866.8 … 311
## 5 5 Copper … 2 (((561012 1148301, 559393.7 11481… 2294
## 6 6 Kodiak 3 (((115112.5 983293, 113051.3 9828… 8126
#plot(pop_region_3338)
pop_mgmt <- pop_region_3338 %>%
group_by(mgmt_area) %>%
summarise(total_pop = sum(total_pop), do_union = F)
plot(pop_mgmt["total_pop"])
rivers_3338 <- read_sf("data/shapefiles/ak_rivers_simp.shp")
st_crs(rivers_3338)
## Coordinate Reference System:
## No EPSG code
## proj4string: "+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"
head(rivers_3338)
## Simple feature collection with 6 features and 3 fields
## geometry type: MULTILINESTRING
## dimension: XY
## bbox: xmin: -622169.8 ymin: 557375.7 xmax: 738652.5 ymax: 2316718
## epsg (SRID): NA
## proj4string: +proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs
## # A tibble: 6 x 4
## StrOrder region n geometry
## <int> <chr> <int> <MULTILINESTRING [m]>
## 1 3 Aleutian I… 316 ((-616509.9 557976.9, -616240.8 557375.7), (-…
## 2 3 Arctic 10869 ((-20959.19 2009003, -20804.35 2009054), (-20…
## 3 3 Bristol Bay 4884 ((-231927.5 794449.9, -232004.6 794730.2), (-…
## 4 3 Chignik 141 ((-302874.8 715936.3, -302829.4 715920.8), (-…
## 5 3 Cook Inlet 3447 ((-11441.06 987666.2, -11772.61 988142.6), (-…
## 6 3 Copper Riv… 1784 ((540428.6 1180052, 539410.1 1176802), (54043…
ggplot() +
geom_sf(data = pop_region_3338, mapping = aes(fill = total_pop)) +
geom_sf(data = pop_3338, mapping = aes(), size = .5) +
geom_sf(data = rivers_3338, mapping = aes(size = StrOrder), color = "Black") +
scale_size(range = c(.01, .2), guide = F) +
theme_bw() +
labs(fill = "Total Population") +
scale_fill_continuous(low = "khaki", high = "firebrick", labels = comma)
pop_3857 <- pop_3338 %>%
st_transform(crs = 3857)
# Define a function to fix the bbox to be in EPSG:3857
# See https://github.com/dkahle/ggmap/issues/160#issuecomment-397055208
ggmap_bbox_to_3857 <- function(map) {
if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
# Extract the bounding box (in lat/lon) from the ggmap to a numeric vector,
# and set the names to what sf::st_bbox expects:
map_bbox <- setNames(unlist(attr(map, "bb")),
c("ymin", "xmin", "ymax", "xmax"))
# Coonvert the bbox to an sf polygon, transform it to 3857,
# and convert back to a bbox (convoluted, but it works)
bbox_3857 <- st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs = 4326)), 3857))
# Overwrite the bbox of the ggmap object with the transformed coordinates
attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
map
}
bbox <- c(-170, 52, -130, 64)
ak_map <- get_stamenmap(bbox, zoom = 4)
## Source : http://tile.stamen.com/terrain/4/0/4.png
## Source : http://tile.stamen.com/terrain/4/1/4.png
## Source : http://tile.stamen.com/terrain/4/2/4.png
## Source : http://tile.stamen.com/terrain/4/0/5.png
## Source : http://tile.stamen.com/terrain/4/1/5.png
## Source : http://tile.stamen.com/terrain/4/2/5.png
ak_map_3857 <- ggmap_bbox_to_3857(ak_map)
ggmap(ak_map_3857) +
geom_sf(data = pop_3857, aes(color = population), inherit.aes = FALSE) +
scale_color_continuous(low = "khaki", high = "firebrick", labels = comma)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
epsg3338 <- leaflet::leafletCRS(
crsClass = "L.Proj.CRS",
code = "EPSG:3338",
proj4def = "+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs",
resolutions = 2^(16:7))
pop_region_4326 <- pop_region_3338 %>%
st_transform(4326)
leaflet(options = leafletOptions(crs = epsg3338)) %>%
addPolygons(data = pop_region_4326, fill = "gray", weight = 1)
pal <- colorNumeric(palette = "Reds", domain = pop_region_4326$total_pop)
m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
addPolygons(data = pop_region_4326,
fillColor = ~pal(total_pop),
weight = 1,
color = "black",
fillOpacity = 1,
label = ~region) %>%
addLegend(position = "bottomleft",
pal = pal,
values = range(pop_region_4326$total_pop),
title = "Total Population")
m
pal <- colorNumeric(palette = "Reds", domain = pop_region_4326$total_pop)
m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
addPolygons(data = pop_region_4326,
fillColor = ~pal(total_pop),
weight = 1,
color = "black",
fillOpacity = 1) %>%
addCircleMarkers(data = pop_4326,
lat = ~lat,
lng = ~lng,
radius = ~log(population/500), # arbitrary scaling
fillColor = "gray",
fillOpacity = 1,
weight = 0.25,
color = "black",
label = ~paste0(pop_4326$city, ", population ", comma(pop_4326$population))) %>%
addLegend(position = "bottomleft",
pal = pal,
values = range(pop_region_4326$total_pop),
title = "Total Population")
m